library(curl)
library(tidyverse)
library(lubridate, warn.conflicts = FALSE)
library(broom.mixed)
library(ggplot2)
library(gridExtra)
library(tidyr)
library(MCMCglmm)
library(lme4)
library(ggplot2)
f <- curl("https://raw.githubusercontent.com/rhottensomers/reesehs-data-replication-assignment/main/doi_10.5061_dryad.kkwh70s31__v4%202/Gao_Cords_2020_IJP_Dataset1_Daily_data.csv") #imports data about zombies
d1 <- read.csv(f, header = TRUE, sep = ",", stringsAsFactors = FALSE) #reads data and creates data frame from it
head(d1) #returns data frame
## Date Group N.sex.active.females N.males.in.group Mating.season
## 1 1/1/06 Gn ND ND NO
## 2 1/2/06 Gn 0 1 NO
## 3 1/3/06 Gn 2 1 NO
## 4 1/4/06 Gn 0 1 NO
## 5 1/5/06 Gn 1 1 NO
## 6 1/6/06 Gn 0 1 NO
## Conception.period Female.Group.Size
## 1 NO 12
## 2 NO 12
## 3 NO 12
## 4 NO 12
## 5 NO 12
## 6 NO 12
str(d1)
## 'data.frame': 17337 obs. of 7 variables:
## $ Date : chr "1/1/06" "1/2/06" "1/3/06" "1/4/06" ...
## $ Group : chr "Gn" "Gn" "Gn" "Gn" ...
## $ N.sex.active.females: chr "ND" "0" "2" "0" ...
## $ N.males.in.group : chr "ND" "1" "1" "1" ...
## $ Mating.season : chr "NO" "NO" "NO" "NO" ...
## $ Conception.period : chr "NO" "NO" "NO" "NO" ...
## $ Female.Group.Size : int 12 12 12 12 12 12 12 12 12 12 ...
names(d1)#mating season and conception period and group as factor, n of sex fs and males in group switch to NA and then convert to factor
## [1] "Date" "Group" "N.sex.active.females"
## [4] "N.males.in.group" "Mating.season" "Conception.period"
## [7] "Female.Group.Size"
d1$N.sex.active.females <- na_if(d1$N.sex.active.females,"ND")
d1$N.males.in.group <- na_if(d1$N.males.in.group,"ND")
d1$Mating.season <- na_if(d1$Mating.season, "ND")
#d1$Group <- as.factor(d1$Group)
d1$Mating.season <- as.factor(d1$Mating.season)
d1$Conception.period <- as.factor(d1$Conception.period)
d1$N.sex.active.females <- as.numeric(d1$N.sex.active.females)
d1$N.males.in.group <- as.numeric(d1$N.males.in.group)
d1$Date <- mdy(d1$Date)
str(d1)
## 'data.frame': 17337 obs. of 7 variables:
## $ Date : Date, format: "2006-01-01" "2006-01-02" ...
## $ Group : chr "Gn" "Gn" "Gn" "Gn" ...
## $ N.sex.active.females: num NA 0 2 0 1 0 0 1 NA NA ...
## $ N.males.in.group : num NA 1 1 1 1 1 1 1 NA NA ...
## $ Mating.season : Factor w/ 2 levels "NO","YES": 1 1 1 1 1 1 1 1 1 1 ...
## $ Conception.period : Factor w/ 2 levels "NO","YES": 1 1 1 1 1 1 1 1 1 1 ...
## $ Female.Group.Size : int 12 12 12 12 12 12 12 12 12 12 ...
view(d1)
mmd <- sum(d1$N.males.in.group >= 2, na.rm = TRUE)
mmd/nrow(d1) * 100 #percent of days with multiple male
## [1] 14.61614
view(d1)
#variables needed for GLMMs
#"test whether larger groups were more likely to contain multiple males and multiple sexually active females per day across group-specific conception periods, we used binomial-family generalized linear mixed models (GLMMs) with a logit- link function and group as a random effect."
#predictor: mean daily f group size
#dependent: proportion of days(during conception period) when a group contain either: more than two males or more than two sexually active females
d1 <- d1 %>% mutate(Mbin = case_when(N.males.in.group < 2 ~ "<2",N.males.in.group >= 2 ~ "2+"))
summary(as.factor(d1$Mbin))
## <2 2+ NA's
## 14185 2534 618
#d1 <- filter(d1, Conception.period == "YES" )
head(d1)
## Date Group N.sex.active.females N.males.in.group Mating.season
## 1 2006-01-01 Gn NA NA NO
## 2 2006-01-02 Gn 0 1 NO
## 3 2006-01-03 Gn 2 1 NO
## 4 2006-01-04 Gn 0 1 NO
## 5 2006-01-05 Gn 1 1 NO
## 6 2006-01-06 Gn 0 1 NO
## Conception.period Female.Group.Size Mbin
## 1 NO 12 <NA>
## 2 NO 12 <2
## 3 NO 12 <2
## 4 NO 12 <2
## 5 NO 12 <2
## 6 NO 12 <2
d1 <- d1 %>% mutate(Fbin = case_when(N.sex.active.females < 2 ~ "<2F",N.sex.active.females >= 2 ~ "2+F"))
head(d1)
## Date Group N.sex.active.females N.males.in.group Mating.season
## 1 2006-01-01 Gn NA NA NO
## 2 2006-01-02 Gn 0 1 NO
## 3 2006-01-03 Gn 2 1 NO
## 4 2006-01-04 Gn 0 1 NO
## 5 2006-01-05 Gn 1 1 NO
## 6 2006-01-06 Gn 0 1 NO
## Conception.period Female.Group.Size Mbin Fbin
## 1 NO 12 <NA> <NA>
## 2 NO 12 <2 <2F
## 3 NO 12 <2 2+F
## 4 NO 12 <2 <2F
## 5 NO 12 <2 <2F
## 6 NO 12 <2 <2F
view(d1)
# a: 2006-06-26- 2006-08-16, Gn
df2 <- d1[c(177:228),]
view(df2)
names(df2)
## [1] "Date" "Group" "N.sex.active.females"
## [4] "N.males.in.group" "Mating.season" "Conception.period"
## [7] "Female.Group.Size" "Mbin" "Fbin"
fdf <- data.frame("Gn", mean(df2$Female.Group.Size), sum(df2$Mbin == "2+")/length(df2$Conception.period), sum(df2$Fbin == "2+F")/length(df2$Conception.period))
names(fdf) <- c("GroupNames", "Mean.Female.Group.Size", "Proportion.of.Two.Plus.Male.Days", "Proportion.of.Two.Plus.SA.Female.Days")
head(fdf)
## GroupNames Mean.Female.Group.Size Proportion.of.Two.Plus.Male.Days
## 1 Gn 12 0.01923077
## Proportion.of.Two.Plus.SA.Female.Days
## 1 0.1538462
calc <- function(f) {
m <- (mean(f$Female.Group.Size))
s1 <- (sum(f$Mbin == "2+", na.rm = TRUE)/length(f$Conception.period))
s2 <- (sum(f$Fbin == "2+F", na.rm = TRUE)/length(f$Conception.period))
n <- f[1,2]
new_row = c(n, m, s1, s2)
return(new_row)
}
# b: 2007-06-09 - 2007-12-15, Gn
b <- d1[c(525:714),]
fdf <- rbind(fdf, calc(b))
# c: 2008-06-18 - 2008-11-01
c <- d1[c(900:1036),]
fdf <- rbind(fdf, calc(c))
# d: 2009-06-13- 2009-10-03
d <- d1[c(1260:1372),]
fdf <- rbind(fdf, calc(d))
# e: 2010-07-03 - 2010-09-20
e <- d1[c(1645:1724),]
fdf <- rbind(fdf, calc(e))
# f: 2012-01-20 - 2012-02-17
f <- d1[c(2211:2239),]
fdf <- rbind(fdf, calc(f))
# g: 2012-07-23 - 2012-11-26
g <- d1[c(2396:2522),]
fdf <- rbind(fdf, calc(g))
# h: 2013-03-01 - 2013-04-05
h <- d1[c(2617:2652),]
fdf <- rbind(fdf, calc(h))
# i: 2013-06-20 - 2013-11-22
i <- d1[c(2728:2883),]
fdf <- rbind(fdf, calc(i))
# j: 2014-04-25 - 2014-11-01
j <- d1[c(3037:3227),]
fdf <- rbind(fdf, calc(j))
# k: 2006-06-10 - 2006-11-11, Gs
k <- d1[c(3448:3602),]
fdf <- rbind(fdf, calc(k))
# l: 2007-06-11 - 2007-09-01
l <- d1[c(3814:3896),]
fdf <- rbind(fdf, calc(l))
# m: 2007-12-11 - 2008-01-08
m <- d1[c(3997:4025),]
fdf <- rbind(fdf, calc(m))
# n: 2008-06-30 - 2008-09-08
n <- d1[c(4199:4269),]
fdf <- rbind(fdf, calc(n))
# o: 2008-10-28 - 2008-11-25, Gsa
o <- d1[c(4319:4347),]
fdf <- rbind(fdf, calc(o))
# p: 2009-06-05 - 2009-07-16
p <- d1[c(4539:4580),]
fdf <- rbind(fdf, calc(p))
# q: 2009-10-11 - 2009-10-29, Gsaa
q <- d1[c(4667:4685),]
fdf <- rbind(fdf, calc(q))
# r: 2010-04-03 - 2010-05-01
r <- d1[c(4841:4869),]
fdf <- rbind(fdf, calc(r))
# s: 2010-06-21 - 2010-08-26
s <- d1[c(4920:4986),]
fdf <- rbind(fdf, calc(s))
# t: 2011-02-10 - 2011-12-13
t <- d1[c(5154:5460),]
fdf <- rbind(fdf, calc(t))
# u: 2012-07-21 - 2012-11-17
u <- d1[c(5681:5800),]
fdf <- rbind(fdf, calc(u))
# v: 2013-05-19 - 2013-09-25
v <- d1[c(5983:6112),]
fdf <- rbind(fdf, calc(v))
# w: 2014-03-06 - 2014-05-27
w <- d1[c(6274:6356),]
fdf <- rbind(fdf, calc(w))
# x: 2014-07-31 - 2014-08-28
x <- d1[c(6421:6449),]
fdf <- rbind(fdf, calc(x))
# y: 2008-10-03 - 2008-11-02, Gsb
y <- d1[c(6575:6605),]
fdf <- rbind(fdf, calc(y))
# z: 2010-03-25 - 2010-04-22
z <- d1[c(7113:7141),]
fdf <- rbind(fdf, calc(z))
# a1: 2010-06-28 - 2010-07-26
a1 <- d1[c(7208:7236),]
fdf <- rbind(fdf, calc(a1))
# b1: 2011-06-09 - 2011-07-07
b1 <- d1[c(7554:7582),]
fdf <- rbind(fdf, calc(b1))
# c1: 2011-12-01 - 2011-12-29
c1 <- d1[c(7729:7757),]
fdf <- rbind(fdf, calc(c1))
# d1: 2012-05-01 - 2012-05-29
done <- d1[c(7881:7909),]
fdf <- rbind(fdf, calc(done))
# e1: 2013-06-08 - 2013-07-06
e1 <- d1[c(8284:8312),]
fdf <- rbind(fdf, calc(e1))
# f1: 2014-02-23 - 2014-03-23
f1 <- d1[c(8544:8572),]
fdf <- rbind(fdf, calc(f1))
# g1: 2010-06-20 - 2010-08-21, Gsc
g1 <- d1[c(9108:9170),]
fdf <- rbind(fdf, calc(g1))
# h1: 2010-10-24 - 2010-11-21
h1 <- d1[c(9234:9262),]
fdf <- rbind(fdf, calc(h1))
# i1: 2011-01-27 - 2011-02-24
i1 <- d1[c(9329:9357),]
fdf <- rbind(fdf, calc(i1))
# j1: 2011-05-11 - 2011-07-15
j1 <- d1[c(9433:9498),]
fdf <- rbind(fdf, calc(j1))
# k1: 2011-12-12 - 2012-01-09
k1 <- d1[c(9648:9676),]
fdf <- rbind(fdf, calc(k1))
# l1: 2012-05-26 - 2012-10-13
l1 <- d1[c(9814:9954),]
fdf <- rbind(fdf, calc(l1))
# m1: 2013-03-01 - 2013-03-29
m1 <- d1[c(10093:10121),]
fdf <- rbind(fdf, calc(m1))
# n1: 2013-05-07 - 2013-08-21
n1 <- d1[c(10160:10266),]
fdf <- rbind(fdf, calc(n1))
# o1: 2014-08-15 - 2014-09-26
o1 <- d1[c(10625:10667),]
fdf <- rbind(fdf, calc(o1))
# p1: 2006-07-14 - 2006-08-15, Twn
p1 <- d1[c(10958:10990),]
fdf <- rbind(fdf, calc(p1))
# q1: 2007-10-05 - 2007-11-02
q1 <- d1[c(11406:11434),]
fdf <- rbind(fdf, calc(q1))
# r1: 2008-07-15 - 2008-09-04
r1 <- d1[c(11690:11741),]
fdf <- rbind(fdf, calc(r1))
# s1: 2009-08-01 - 2009-10-23
s1 <- d1[c(12072:12155),]
fdf <- rbind(fdf, calc(s1))
# t1: 2010-07-10 - 2010-09-12
t1 <- d1[c(12415:12479),]
fdf <- rbind(fdf, calc(t1))
# u1: 2011-04-26 - 2011-06-23
u1 <- d1[c(12705:12763),]
fdf <- rbind(fdf, calc(u1))
# v1: 2011-08-17 - 2011-09-14
v1 <- d1[c(12818:12846),]
fdf <- rbind(fdf, calc(v1))
# w1: 2012-08-04 - 2012-11-07
w1 <- d1[c(13171:13266),]
fdf <- rbind(fdf, calc(w1))
# x1: 2013-03-03 - 2013-03-31
x1 <- d1[c(13382:13410),]
fdf <- rbind(fdf, calc(x1))
# y1: 2013-07-16 - 2013-08-14
y1 <- d1[c(13517:13546),]
fdf <- rbind(fdf, calc(y1))
# z1: 2014-04-03 - 2014-05-01
z1 <- d1[c(13778:13806),]
fdf <- rbind(fdf, calc(z1))
# a2: 2014-06-09 - 2014-10-21
a2 <- d1[c(13845:13979),]
fdf <- rbind(fdf, calc(a2))
# b2: 2006-06-10 - 2006-08-24, Tws
b2 <- d1[c(14211:14286),]
fdf <- rbind(fdf, calc(b2))
# c2: 2007-08-03 - 2007-11-04
c2 <- d1[c(14630:14723),]
fdf <- rbind(fdf, calc(c2))
# d2: 2008-07-10 - 2008-10-15
dtwo <- d1[c(14972:15069),]
fdf <- rbind(fdf, calc(dtwo))
# e2: 2008-12-16 - 2009-02-15
e2 <- d1[c(15131:15192),]
fdf <- rbind(fdf, calc(e2))
# f2: 2010-01-04 - 2010-02-01
f2 <- d1[c(15515:15543),]
fdf <- rbind(fdf, calc(f2))
# g2: 2010-04-26 - 2010-09-15
g2 <- d1[c(15627:15769),]
fdf <- rbind(fdf, calc(g2))
# h2: 2010-12-27 - 2011-02-15
h2 <- d1[c(15872:15922),]
fdf <- rbind(fdf, calc(h2))
# i2: 2011-05-05 - 2011-12-11
i2 <- d1[c(16001:16221),]
fdf <- rbind(fdf, calc(i2))
# j2: 2012-07-28 - 2012-11-23
j2 <- d1[c(16451:16569),]
fdf <- rbind(fdf, calc(j2))
# k2: 2013-05-30 - 2013-10-20
k2 <- d1[c(16757:16900),]
fdf <- rbind(fdf, calc(k2))
# l2: 2014-04-06 - 2014-05-04
l2 <- d1[c(17068:17096),]
fdf <- rbind(fdf, calc(l2))
# m2: 2014-06-25 - 2014-12-31
m2 <- d1[c(17148:17337),]
fdf <- rbind(fdf, calc(m2))
head(fdf)
## GroupNames Mean.Female.Group.Size Proportion.of.Two.Plus.Male.Days
## 1 Gn 12 0.0192307692307692
## 2 Gn 12 0.0421052631578947
## 3 Gn 12.7226277372263 0.0145985401459854
## 4 Gn 13 0.0176991150442478
## 5 Gn 16 0.0375
## 6 Gn 14 0.0689655172413793
## Proportion.of.Two.Plus.SA.Female.Days
## 1 0.153846153846154
## 2 0.110526315789474
## 3 0.0802919708029197
## 4 0.176991150442478
## 5 0.2
## 6 0.586206896551724
fdf$GroupNames <- as.factor(fdf$GroupNames )
names(fdf)
## [1] "GroupNames"
## [2] "Mean.Female.Group.Size"
## [3] "Proportion.of.Two.Plus.Male.Days"
## [4] "Proportion.of.Two.Plus.SA.Female.Days"
fdf$Mean.Female.Group.Size <- as.numeric(fdf$Mean.Female.Group.Size)
fdf$Proportion.of.Two.Plus.Male.Days <- as.numeric(fdf$Proportion.of.Two.Plus.Male.Days)
fdf$Proportion.of.Two.Plus.SA.Female.Days <- as.numeric(fdf$Proportion.of.Two.Plus.SA.Female.Days)
view(fdf)
names(fdf)
## [1] "GroupNames"
## [2] "Mean.Female.Group.Size"
## [3] "Proportion.of.Two.Plus.Male.Days"
## [4] "Proportion.of.Two.Plus.SA.Female.Days"
str(fdf)
## 'data.frame': 65 obs. of 4 variables:
## $ GroupNames : Factor w/ 8 levels "Gn","Gs","Gsa",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Mean.Female.Group.Size : num 12 12 12.7 13 16 ...
## $ Proportion.of.Two.Plus.Male.Days : num 0.0192 0.0421 0.0146 0.0177 0.0375 ...
## $ Proportion.of.Two.Plus.SA.Female.Days: num 0.1538 0.1105 0.0803 0.177 0.2 ...
#Gs
dfgs <- filter(d1, Group == "Gs" )
#Gn
dfgn <- filter(d1, Group == "Gn" )
#Gsaa
dfgsaa <- filter(d1, Group == "Gsaa" )
#Gsa
dfgsa <- filter(d1, Group == "Gsa" )
#Gsc
dfgsc <- filter(d1, Group == "Gsc" )
#Twn
dftwn <- filter(d1, Group == "Twn" )
#Gsb
dfgsb <- filter(d1, Group == "Gsb" )
#Tws
dftws <- filter(d1, Group == "Tws" )
tdf <- data.frame("Tws",sum(dftws$Female.Group.Size)/length(dftws$Female.Group.Size), length(dftws$Group), sum(dftws$Mbin == "2+", na.rm = TRUE)/length(dftws$Mbin), sum(dftws$Conception.period == "YES"))
names(tdf) <- c("Group", "Mean female group size", "Total observation days", "Proportion of multimale days", "Total group-specific conception period days")
dftws <- filter(dftws, Conception.period == "YES" )
tdf$`Proportion of conception period days with.multiple males` <- sum(dftws$Mbin == "2+", na.rm = T)/length(dftws$Mbin)
table <- function(f){
gn <- f[1,2]
mfgs <- sum(f$Female.Group.Size)/length(f$Female.Group.Size)
od <- length(f$Group)
pmmd <- sum(f$Mbin == "2+", na.rm = TRUE)/length(f$Group)
cpd <- sum(f$Conception.period == "YES" )
f <- filter(f, Conception.period == "YES" )
mmcpd <- sum(f$Mbin == "2+", na.rm = TRUE)/length(f$Mbin)
new_dfrow = c(gn, mfgs, od, pmmd, cpd, mmcpd)
return(new_dfrow)
}
tdf <- rbind(tdf, table(dfgs))
tdf <- rbind(tdf, table(dfgn))
tdf <- rbind(tdf, table(dfgsaa))
tdf <- rbind(tdf, table(dfgsa))
tdf <- rbind(tdf, table(dfgsc))
tdf <- rbind(tdf, table(dftwn))
tdf <- rbind(tdf, table(dfgsb))
head(tdf, n = 8)
## Group Mean female group size Total observation days
## 1 Tws 17.2303011864922 3287
## 2 Gs 16.8300198807157 1006
## 3 Gn 13.6632187404928 3287
## 4 Gsaa 10.7578616352201 1908
## 5 Gsa 9 373
## 6 Gsc 8.35849056603774 1908
## 7 Twn 7.79738363249163 3287
## 8 Gsb 3.40070144673389 2281
## Proportion of multimale days Total group-specific conception period days
## 1 0.330088226346212 1256
## 2 0.175944333996024 338
## 3 0.122908427137207 1111
## 4 0.0849056603773585 784
## 5 0.0723860589812333 71
## 6 0.0178197064989518 536
## 7 0.194706419227259 670
## 8 0.00219202104340202 234
## Proportion of conception period days with.multiple males
## 1 0.493630573248408
## 2 0.233727810650888
## 3 0.183618361836184
## 4 0.151785714285714
## 5 0.0704225352112676
## 6 0.0317164179104478
## 7 0.12089552238806
## 8 0
str(tdf)
## 'data.frame': 8 obs. of 6 variables:
## $ Group : chr "Tws" "Gs" "Gn" "Gsaa" ...
## $ Mean female group size : chr "17.2303011864922" "16.8300198807157" "13.6632187404928" "10.7578616352201" ...
## $ Total observation days : chr "3287" "1006" "3287" "1908" ...
## $ Proportion of multimale days : chr "0.330088226346212" "0.175944333996024" "0.122908427137207" "0.0849056603773585" ...
## $ Total group-specific conception period days : chr "1256" "338" "1111" "784" ...
## $ Proportion of conception period days with.multiple males: chr "0.493630573248408" "0.233727810650888" "0.183618361836184" "0.151785714285714" ...
tdf$`Mean female group size` <- as.numeric(tdf$`Mean female group size`)
tdf$`Proportion of multimale days` <- as.numeric(tdf$`Proportion of multimale days`)
tdf$`Total group-specific conception period days` <- as.numeric(tdf$`Total group-specific conception period days`)
tdf$`Proportion of conception period days with.multiple males` <- as.numeric(tdf$`Proportion of conception period days with.multiple males`)
tdf$`Mean female group size` <- round(tdf$`Mean female group size`, digits = 1)
tdf$`Proportion of multimale days` <- round(tdf$`Proportion of multimale days`, digits = 3)
tdf$`Proportion of conception period days with.multiple males` <- round(tdf$`Proportion of conception period days with.multiple males`, digits = 2)
knitr::kable(tdf)
| Group | Mean female group size | Total observation days | Proportion of multimale days | Total group-specific conception period days | Proportion of conception period days with.multiple males |
|---|---|---|---|---|---|
| Tws | 17.2 | 3287 | 0.330 | 1256 | 0.49 |
| Gs | 16.8 | 1006 | 0.176 | 338 | 0.23 |
| Gn | 13.7 | 3287 | 0.123 | 1111 | 0.18 |
| Gsaa | 10.8 | 1908 | 0.085 | 784 | 0.15 |
| Gsa | 9.0 | 373 | 0.072 | 71 | 0.07 |
| Gsc | 8.4 | 1908 | 0.018 | 536 | 0.03 |
| Twn | 7.8 | 3287 | 0.195 | 670 | 0.12 |
| Gsb | 3.4 | 2281 | 0.002 | 234 | 0.00 |
For this data replication I decided to replicate two GLMMs that described whether a particular group was more likely to include multiple(more than two) males or multiple sexually active females during the already defined 65 conception periods. This required two GLMMs: one that looked at multiple males and group size, and another that looked at multiple sexual active females and group size. The GLMMs themselves as states within the statistics section of the article were binomial family and thus used a logit-link function and group as a random effect. The predictor values being mean female group size during each conception period and proportion of days during conception periods where there were multiple males or sexual active females However, it was not stated whether they used intercept or slope form for their GLMMs. So I decided to run both to see if one fit the data better, and the intercept form resulted with a better AIC(52.2 vs 56.2). I will first begin with the multiple male GLMM below:
In the results section of the article it was stated that their GLMM resulted in these values: β = 0.021, 95% CI [0.007, 0.034], P = 0.015. In this section I will run the GLMM and attempt to receive the same or similar numbers.
mglmm <-glmer( Proportion.of.Two.Plus.Male.Days ~ Mean.Female.Group.Size + (1 | GroupNames), family = binomial, data = fdf) #intercept form of multi-male GLMM with binomial link and fdf data frame of conception periods
summary(mglmm) #returns glmm
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Proportion.of.Two.Plus.Male.Days ~ Mean.Female.Group.Size + (1 |
## GroupNames)
## Data: fdf
##
## AIC BIC logLik deviance df.resid
## 52.2 58.7 -23.1 46.2 62
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5171 -0.3147 -0.2199 0.1676 2.8253
##
## Random effects:
## Groups Name Variance Std.Dev.
## GroupNames (Intercept) 0 0
## Number of obs: 65, groups: GroupNames, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.38492 1.10757 -3.056 0.00224 **
## Mean.Female.Group.Size 0.11846 0.07985 1.483 0.13795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Mn.Fml.Gr.S -0.938
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
null <- glmer(Proportion.of.Two.Plus.Male.Days ~ 1 + (1|GroupNames), family = binomial(link = "logit"), data = fdf) #reduced or null version of the GLMM used for liklihood ratio test
anova(null, mglmm, test = "Chisq") #liklihood ratio test and returns p value
## Data: fdf
## Models:
## null: Proportion.of.Two.Plus.Male.Days ~ 1 + (1 | GroupNames)
## mglmm: Proportion.of.Two.Plus.Male.Days ~ Mean.Female.Group.Size + (1 | GroupNames)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## null 2 52.491 56.840 -24.246 48.491
## mglmm 3 52.182 58.706 -23.091 46.182 2.309 1 0.1286
exp(confint(mglmm, parm = "Mean.Female.Group.Size")) #returns 95% confidence intervals
## 2.5 % 97.5 %
## Mean.Female.Group.Size 0.9665444 1.331562
Discussion, discussion
smglmm <- glmer( Proportion.of.Two.Plus.Male.Days ~ Mean.Female.Group.Size + (1 + Mean.Female.Group.Size | GroupNames), family = binomial(link = "logit"), data = fdf) #slope form of multi-male GLMM
summary(smglmm) #returns glmm
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Proportion.of.Two.Plus.Male.Days ~ Mean.Female.Group.Size + (1 +
## Mean.Female.Group.Size | GroupNames)
## Data: fdf
##
## AIC BIC logLik deviance df.resid
## 56.2 67.1 -23.1 46.2 60
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5171 -0.3147 -0.2199 0.1676 2.8253
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## GroupNames (Intercept) 4.498e-10 2.121e-05
## Mean.Female.Group.Size 1.062e-12 1.031e-06 -1.00
## Number of obs: 65, groups: GroupNames, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.38491 1.10758 -3.056 0.00224 **
## Mean.Female.Group.Size 0.11846 0.07985 1.483 0.13796
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Mn.Fml.Gr.S -0.938
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
AIC values etc
mlmer <- lmer(data = fdf, Proportion.of.Two.Plus.Male.Days ~ Mean.Female.Group.Size + (1|GroupNames), REML = FALSE) #linear mixed model version of mutli-male GLMM
summary(mlmer) #returns multi-male lmm
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Proportion.of.Two.Plus.Male.Days ~ Mean.Female.Group.Size + (1 |
## GroupNames)
## Data: fdf
##
## AIC BIC logLik deviance df.resid
## 10.2 18.9 -1.1 2.2 61
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.29632 -0.61023 -0.17619 0.03434 3.01982
##
## Random effects:
## Groups Name Variance Std.Dev.
## GroupNames (Intercept) 0.00000 0.0000
## Residual 0.06059 0.2461
## Number of obs: 65, groups: GroupNames, 8
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.069055 0.075004 -0.921
## Mean.Female.Group.Size 0.021780 0.006258 3.480
##
## Correlation of Fixed Effects:
## (Intr)
## Mn.Fml.Gr.S -0.913
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
coefficients(mlmer)
## $GroupNames
## (Intercept) Mean.Female.Group.Size
## Gn -0.06905499 0.02178015
## Gs -0.06905499 0.02178015
## Gsa -0.06905499 0.02178015
## Gsaa -0.06905499 0.02178015
## Gsb -0.06905499 0.02178015
## Gsc -0.06905499 0.02178015
## Twn -0.06905499 0.02178015
## Tws -0.06905499 0.02178015
##
## attr(,"class")
## [1] "coef.mer"
nmlmer <- lmer(data = fdf, Proportion.of.Two.Plus.Male.Days ~ 1 + (1|GroupNames), REML = FALSE) #reduced form of linear mixed model
anova(nmlmer, mlmer, test = "Chisq") #likelihood ratio test
## Data: fdf
## Models:
## nmlmer: Proportion.of.Two.Plus.Male.Days ~ 1 + (1 | GroupNames)
## mlmer: Proportion.of.Two.Plus.Male.Days ~ Mean.Female.Group.Size + (1 | GroupNames)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nmlmer 3 14.188 20.711 -4.0938 8.1877
## mlmer 4 10.221 18.919 -1.1105 2.2210 5.9667 1 0.01458 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(mlmer, level = .95)
## 2.5 % 97.5 %
## .sig01 0.000000000 0.1186721
## .sigma 0.208732922 0.2953741
## (Intercept) -0.218256767 0.1107460
## Mean.Female.Group.Size 0.006038569 0.0342281
Also in the results section the multi-sexually active female GLMM results were given: β = 0.010, 95% CI [0.002, 0.018], P = 0.017. In this next section I will again attempt to recieve similar numbers from my own GLMMs.
fglmm <- glmer( Proportion.of.Two.Plus.SA.Female.Days ~ Mean.Female.Group.Size + (1 | GroupNames), family =binomial(link = "logit"), data = fdf) #
summary(fglmm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Proportion.of.Two.Plus.SA.Female.Days ~ Mean.Female.Group.Size +
## (1 | GroupNames)
## Data: fdf
##
## AIC BIC logLik deviance df.resid
## 35.9 42.4 -14.9 29.9 62
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.29137 -0.04803 0.19523 0.60640 2.84704
##
## Random effects:
## Groups Name Variance Std.Dev.
## GroupNames (Intercept) 0 0
## Number of obs: 65, groups: GroupNames, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.25538 1.35053 -2.410 0.0159 *
## Mean.Female.Group.Size 0.04642 0.10461 0.444 0.6573
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Mn.Fml.Gr.S -0.924
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
nflgmm <- glmer( Proportion.of.Two.Plus.SA.Female.Days ~ 1 + (1|GroupNames), family =binomial(link = "logit"), data = fdf)#reduced version
anova(nflgmm, fglmm, test = "Chisq")#likelihood ratio test fr
## Data: fdf
## Models:
## nflgmm: Proportion.of.Two.Plus.SA.Female.Days ~ 1 + (1 | GroupNames)
## fglmm: Proportion.of.Two.Plus.SA.Female.Days ~ Mean.Female.Group.Size + (1 | GroupNames)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nflgmm 2 34.053 38.402 -15.027 30.053
## fglmm 3 35.857 42.380 -14.929 29.857 0.1961 1 0.6579
exp(confint(fglmm, parm = "Mean.Female.Group.Size")) # 95% confidence interval
## 2.5 % 97.5 %
## Mean.Female.Group.Size 0.847016 1.295848
wow, this one is weird too!
glmer( Proportion.of.Two.Plus.SA.Female.Days ~ Mean.Female.Group.Size + (1 + Mean.Female.Group.Size | GroupNames), family =binomial(link = "logit"), data = fdf)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Proportion.of.Two.Plus.SA.Female.Days ~ Mean.Female.Group.Size +
## (1 + Mean.Female.Group.Size | GroupNames)
## Data: fdf
## AIC BIC logLik deviance df.resid
## 39.8572 50.7292 -14.9286 29.8572 60
## Random effects:
## Groups Name Std.Dev. Corr
## GroupNames (Intercept) 1.066e-05
## Mean.Female.Group.Size 6.346e-07 -1.00
## Number of obs: 65, groups: GroupNames, 8
## Fixed Effects:
## (Intercept) Mean.Female.Group.Size
## -3.25541 0.04642
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
AIC yadaydad
flmm <- lmer(data = fdf, Proportion.of.Two.Plus.SA.Female.Days ~ Mean.Female.Group.Size + (1|GroupNames), REML = FALSE)
summary(flmm)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Proportion.of.Two.Plus.SA.Female.Days ~ Mean.Female.Group.Size +
## (1 | GroupNames)
## Data: fdf
##
## AIC BIC logLik deviance df.resid
## -52.6 -43.9 30.3 -60.6 61
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4661 -0.5573 -0.2981 0.3277 3.6237
##
## Random effects:
## Groups Name Variance Std.Dev.
## GroupNames (Intercept) 0.00000 0.0000
## Residual 0.02303 0.1518
## Number of obs: 65, groups: GroupNames, 8
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.042139 0.046247 0.911
## Mean.Female.Group.Size 0.010610 0.003859 2.750
##
## Correlation of Fixed Effects:
## (Intr)
## Mn.Fml.Gr.S -0.913
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
nflmm <- lmer(data = fdf, Proportion.of.Two.Plus.SA.Female.Days ~ 1 + (1|GroupNames), REML = FALSE)
anova(nflmm, flmm, test = "Chisq")
## Data: fdf
## Models:
## nflmm: Proportion.of.Two.Plus.SA.Female.Days ~ 1 + (1 | GroupNames)
## flmm: Proportion.of.Two.Plus.SA.Female.Days ~ Mean.Female.Group.Size + (1 | GroupNames)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nflmm 3 -49.232 -42.709 27.616 -55.232
## flmm 4 -52.640 -43.942 30.320 -60.640 5.4076 1 0.02005 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(flmm,level = .95) # 95% confidence interval
## 2.5 % 97.5 %
## .sig01 0.000000000 0.06587508
## .sigma 0.128994523 0.18212604
## (Intercept) -0.049850626 0.14095299
## Mean.Female.Group.Size 0.002173912 0.01828766
fdf$prediction1 <- predict(mglmm, type = "response", newdata = fdf) #adds prediction line using the mutli-male GLMM
p1 <- ggplot(data = fdf, aes(x = Mean.Female.Group.Size, y = Proportion.of.Two.Plus.Male.Days, shape = GroupNames)) + geom_point() + scale_shape_manual( values = c(2,1,4,3,7,5,6,0)) + geom_abline(slope = .02178, intercept = -0.06905, color = "red") + ylab("Proportion of Multimale Days") + xlab("Female Group Size") + geom_line(data = fdf, aes(x = Mean.Female.Group.Size, y = prediction1), color = "blue") + scale_color_manual(name = "Models", values = c("Linear Mixed Model" = "red", "GLMM" = "blue")) +theme(legend.direction = "vertical", legend.box = "vertical") #plots conception period data along with the linear mixed model and the GLMM lines
fdf$prediction2 <- predict(fglmm, type = "response", newdata = fdf)#adds prediction line for mutli-sa female GLMM
p2 <- ggplot(data = fdf, aes(x = Mean.Female.Group.Size, y = Proportion.of.Two.Plus.SA.Female.Days, shape = GroupNames )) + geom_point() + ylim(0, 1) + scale_shape_manual( values = c(2,1,4,3,7,5,6,0)) + geom_abline(slope = 0.01061 , intercept = 0.04214, color = "red" ) + ylab("Proportion of Days with ≥2 Sexual Active Females") + xlab("Female Group Size") + geom_line(data = fdf, aes(x = Mean.Female.Group.Size, y = prediction2), color = "blue")
p1
p2